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Abstract 

An added-mass partitioned (AMP) algorithm is described for solving fluid-structure interaction (FSI) 
problems coupling incompressible flows with thin elastic structures undergoing hnite deformations. The new 
AMP scheme is fully second-order accurate and stable, without sub-time-step iterations, even for very light 
structures when added-mass effects are strong. The fluid, governed by the incompressible Navier-Stokes 
equations, is solved in velocity-pressure form using a fractional-step method; large deformations are treated 
with a mixed Eulerian-Lagrangian approach on deforming composite grids. The motion of the thin structure 
is governed by a generalized Euler-Bernoulli beam model, and these equations are solved in a Lagrangian 
frame using two approaches, one based on finite differences and the other on finite elements. The key AMP 
interface condition is a generalized Robin (mixed) condition on the fluid pressure. This condition, which is 
derived at a continuous level, has no adjustable parameters and is applied at the discrete level to couple the 
partitioned domain solvers. Special treatment of the AMP condition is required to couple the finite-element 
beam solver with the finite-difference-based fluid solver, and two coupling approaches are described. A 
normal-mode stability analysis is performed for a linearized model problem involving a beam separating two 
fluid domains, and it is shown that the AMP scheme is stable independent of the ratio of the mass of the fluid 
to that of the structure. A traditional partitioned (TP) scheme using a Dirichlet-Neumann coupling for the 
same model problem is shown to be unconditionally unstable if the added mass of the fluid is too large. A 
series of benchmark problems of increasing complexity are considered to illustrate the behavior of the AMP 
algorithm, and to compare the behavior with that of the TP scheme. The results of all these benchmark 
problems verify the stability and accuracy of the AMP scheme. Results for one benchmark problem modeling 
blood flow in a deforming artery are also compared with corresponding results available in the literature. 
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1. Introduction 


Fluid-structure interaction (FSI) problems that describe the motion of an incompressible fluid coupled to 
a thin walled structure (beam or shell) arise in many applications such as those in structural engineering and 
biomedicine. Such FSI problems are often modeled mathematically by suitable partial differential equations 
for the fluids and structures in their respective domains, together with matching conditions at the boundaries 
of the domains where the solutions of the equations interact. Numerical algorithms used to solve these FSI 
problems can be classified into two main categories. Algorithms belonging to one category, called monolithic 
schemes, treat the equations for the fluids and structures along with interface and boundary conditions 
as a large system of evolution equations, and then advance the solutions together. The other category of 
algorithms are partitioned schemes (also known as modular or sequential schemes), and these algorithms 
employ separate solvers for the fluids and structures which are coupled at the interface. Sub-iterations 
are often performed at each time step of partitioned algorithms for stability. Even though many existing 
partitioned schemes suffer from moderate to serious stability issues in certain problem regimes, they are often 
preferred since they can make use of existing solvers and can be more efficient than monolithic schemes. 

The traditional partitioned algorithm for beams (or shells) uses the velocity and/or acceleration of the 
solid as a boundary condition on the fluid. The force of the fluid is accounted for through a body forcing on 
the beam. It has been found that partitioned schemes may be unstable, or require multiple sub-iterations 
per time step, when the density of the structure is similar to or lighter than that of the fluid [1, 2]. These 
instabilities are attributed to the added-mass effect whereby the force required to accelerate a structure 
immersed in a fluid must also account for accelerating the surrounding fluid. The added-mass effect has 
been found to be especially problematic in many biological flows such as haemodynamics since the density 
of the fluid (blood) is similar to that of the adjacent structure (arterial walls) [3]. 

In our previous companion papers [4, 5], we developed stable added-mass partitioned (AMP) algorithms 
for linearized FSI problems involving incompressible Stokes fluids coupled to elastic bulk solids and to beams. 
The key ingredient of our approach is the use of non-traditional Robin interface conditions. These condition 
are derived at a continuous level, have no adjustable parameters, and are amenable for incorporation in 
high-order accurate schemes. Using mode analysis, the AMP approach was shown to be stable without sub¬ 
iterations per time step, and the numerical results demonstrated second-order accuracy in the max-norm. 
In the present paper, we focus on the case of thin structures and our primary purpose is to extend the AMP 
scheme in [5] to the case of finite amplitude motions in more general geometries. In particular, the AMP 
interface conditions are derived for beams of finite thickness, for beams with fluid on two sides, and for beams 
with a free end immersed in the fluid. Finite amplitude motions of the beam results in finite deformations of 
the surrounding fluid domain. Our numerical approach for the solution of the fluid equations in an evolving 
fluid domain is based on the use of deforming composite grids (DCG). The DCG approach for FSI problems 
was described first in [6] inviscid compressible flow coupled to a linearly elastic solid, and later in [7] for 
compressible flow coupled to nonlinear hyperelastic solids. The approach is extended here for incompressible 
flow coupled to beams. Both finite-difference and finite-element approximations to the equations governing 
the motion of the beam are considered, and issues concerning the interface coupling of a finite-element based 
beam solver to a finite-difference based fluid solver are discussed. The stability analysis in [5] is also extended 
to treat beams with fluid on two sides. Finally, several benchmark FSI problems are presented to illustrate 
the stability and accuracy of the present AMP algorithm. 

The investigation of FSI problems and the development of numerical approximations for their solution 
are very active areas of research, see for example [8-10] and the references therein. For the coupling of incom¬ 
pressible flows and bulk solids, the work in [3, 4, 10-12] describes recent developments. The development of 
partitioned schemes for the case of incompressible fluids coupled to thin structures, the focus of this paper, 
is also an active area. The first stable partitioned scheme for incompressible flow coupled to thin structures 
was the “kinematically coupled scheme” of Guidoboni et al. [13] which uses an operator splitting of the 
kinematic interface condition (matching the fluid and beam velocities). This scheme was further advanced 
in [14, 15]. Lukacova-Medvid’ova et al. [16] developed a partitioned scheme based on the kinematically 
coupled scheme that uses Strang splitting. Fernandez and collaborators have also developed stable parti¬ 
tioned schemes for incompressible flow coupled to thin structures that are based on a time-splitting of the 
kinematic boundary condition [17-22]. However, it appears to be difficult to achieve higher than first-order 
accuracy with the time-splitting approach. The AMP approximation, in contrast, is not based on a time 
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splitting and is amenable to second- or even higher-order accuracy. Second-order accuracy in the max-norm 
was demonstrated in [5] for linearized problems, while in this paper second-order accuracy is shown for the 
full nonlinear problem with deforming domains. 

The remainder of the paper is organized as follows. In Section 2 we describe the governing equations. The 
AMP interface conditions are derived in Section 3. A second-order accurate predictor-corrector algorithm 
based on these conditions is described in Section 4. A specific choice for the beam model is given in Section 5. 
In Section 6, the stability of the AMP scheme is shown for a linearized FSI problem involving a beam with 
fluid on two sides. The numerical approach used for moving domains based on deforming composite grids, our 
numerical approaches for the beam solver, and the treatment of the AMP interface conditions are described 
in Section 7. Numerical results presented in Section 8 carefully demonstrate the stability and accuracy of 
the AMP algorithm. Conclusions are provided in Section 9. 

2. Governing equations 

We consider the fluid-structure coupling of an incompressible fluid and a thin deformable structure. The 
fluid occupies the domain x S r2(t) while the structure lies in the domain x G 12(f), where x is position and 
t is time. The coupling of the fluid and structure occurs along the interface r(f), see Figure 1. It is assumed 
that the fluid is governed by the incompressible Navier-Stokes equations, which in an Eulerian frame are 
given by 

^-t (v • V)v = • (T, xen(f), (1) 

V-v = 0, xen(f), (2) 

where p is the (constant) fluid density and v = v(x, f) is the fluid velocity. The fluid stress tensor, cr = <t(x, f), 
is given by 

(T = -pi + T, T = p [Vv -I- (Vv)"^] , (3) 

where p = p(x, f) is the pressure, I is the identity tensor, r is the viscous stress tensor, and p is the (constant) 
fluid viscosity. For future reference, the components of a vector such as v will be denoted hy Vm, m = 1, 2, 3 
(i.e. V = [z;i, 1 ) 2 ,wa]^), while components of a tensor such as a will be denoted by amn, rn,n = 1,2,3. 
The velocity-divergence form of the equations given by (1) and (2) require appropriate initial and boundary 
conditions, as well as conditions on r(f) where the behavior of the fluid is coupled to that of the solid (as 
discussed below). 



Figure 1: Geometry of a three-dimensional beam immersed in a fluid. The solid domain is Cl{t). The fluid domain is 
f2(f). The interface between the fluid and solid is r(f). The beam reference-line is xo(s,f). A point on the surface of 
the beam is denoted by ^b{6,s,t) where 6 € P(s) is the parameterization of the cross-sectional surface curve. The 
fluid tractions on the surface points X6(0,s,t), for 0 € P(s), contribute to the force on the beam at xo(s,t). 
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An elliptic equation for the fluid pressure can be derived from (l)-(3) as 

Ap= —pVv : (Vv)^, X e Q{t). 

In the velocity-pressure form of the incompressible Navier-Stokes equations, the Poisson equation (4) is used 
in place of (2). This form requires an additional boundary condition, and a natural choice is V • v = 0 for 
X S dn{t), see [23] for example. 

The structure is assumed to be thin and governed by a beam model, which describes the evolution of the 
displacement, u = u(s, t), of a reference centerline curve of the beam as a function of a Lagrangian coordinate 
s G TZ and time t. (Overbars are used throughout to denote quantities belonging to the structure.) The 
beam model has the form 

=L(u,v)-bf(s,t), sgH, (5) 

where p is the (constant) beam density, A = A(s) is the cross-sectional area of the beam, L is the internal 
restoring force per unit length which depends on u and on the velocity, v = 9u/(9t, of the reference curve (in 
general), and f is the external force per unit length on the beam from the fluid (defined below). The effects 
of torsion are neglected in the present beam model. In physical space, the position Xq of tfi® reference curve 
is given by 

Xo(sA) = Xo(s, 0)-b u(s,t), sGTZ, (6) 

where Xo(s, 0) is the initial position of the reference curve of the beam. The position x/, of the surface of the 
beam is given by 

xi,(d,s,t) = xo(s,t) + z(0,s,t), OgP, sgP, (7) 

where z is distance from the reference curve to the perimeter V = V{s) of the cross section of the beam at 
s, parameterized by arclength 6 along the perimeter. Note that z is defined in terms of the beam variables; 
its form for an Euler-Bernoulli beam is given in Section 5. The surface of the beam given by (7) defines the 
position of the fluid-structure interface r(t) (see Figure 1). 

At the fluid-structure interface, the velocity of the fluid at a point on the surface of the beam matches 
the corresponding velocity of the beam surface, 

v(xh(0, s,t),t) = Vb(0, s,t), OgP, sgP, (8) 

where 

d d 

Vb{e,s,t) = —Xbi9,s,t) = -v{s,t) + \v{e,s,t), w{e,s,t) =—z{e,s,t). (9) 

Formulas for the finite-thickness corrections, z and w, are given in Section 5 for the case of an Fuler-Bernoulli 
beam model. The dynamic matching condition (balance of forces) defines the force per unit length on the 
beam, f, in (5) in terms of the fluid force on the surface of the beam, 

f(s,t) = — j {crn){6,s,t) d9, sgP, (10) 

Jp 

where 

(crn)(0,s,t) = (T{xb{9,s,t),t)n{xb{9,s,t),t). (11) 

Here, n denotes the unit normal to the beam surface (outward pointing from the fluid domain). 

We note that for the case of a beam in two dimensions as shown in Figure 2, the surface position in (7) 
reduces to 

^b,±is,t) =Xois,t)+ z±{s,t), sGU, (12) 

where the plus and minus signs denote the upper and lower surfaces of the beam, respectively. The velocity 
on the surfaces of the beam are now given by 

d _ _ _ _ d _ 

Vh,±(sA) = ^Xb,±(s,t) = v(s,t)-b w±(s,t), W±(s,t) = —z±(s,t), (13) 
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Figure 2: Geometry of a two-dimensional beam immersed in a fluid. The beam reference surface is xo(s, t). The top 
and bottom surfaces of the beam are the x+(s,f) and x_(s,t), respectively. 


and if fluid regions exist on both surfaces of the beam, then the the external force is given by 


f(s,t)= (crn)+(s,t)-f (crn)_(s,t) 


s G 7^, 


(14) 


where 


(crn)±(s,t) = cr(xt,_±(s,t),f)n(xb_±(s,t),f). (15) 

The formulas (12)-(15) also hold for the case of a thin plate or shell provided the dependent variables, xo(s, t), 
v(s,t), etc., are taken to depend on two independent parameters (si,S 2 ) = s G that parameterize the 
reference surface of the shell. 


3. Added-Mass Partitioning 

In the traditional approach to partitioned fluid-structure interaction, the kinematic condition (8) is 
satisfied by imposing the motion of the solid as an interface condition on the fluid. This inherently one-sided 
approach is one root cause leading to added-mass instabilities for partitioned fluid-structure interaction. 
Various avenues to more symmetric imposition of (8) have been taken with various degrees of success. Here 
we pursue the Added Mass Partitioned (AMP) approach to coupling. The core AMP algorithm for coupling a 
linear incompressible fluid with a thin shell or beam was described in [5] . In that work a simplified model was 
considered where linearized equations of motion were used to describe infinitesimal amplitude perturbations 
from the specihc geometry of fluid in a box with a beam on one of the faces of the box. This manuscript 
considers the more general case of the nonlinear Navier-Stokes equations with finite amplitude deformations 
applied to more general geometries. 

3.1. Continuous AMP Interface Condition 

The AMP interface condition can be derived at a continuous level by first taking the time derivative of 
the kinematic interface condition (8) and using the definition of the velocity on the surface of the beam (9) 
to give 


D d d 

—v(xb(6»,s,t),t) = ^v(s,t) -b —w{e,s,t) 

d 

= ^{s,t) + ^^Md,s,t), (16) 

where D/Dt = 9* -I- v • V is the material derivative. The fluid acceleration, Dv/Dt^ is given by the fluid 
momentum equation (1), while the acceleration of the structure, is given by the beam equation (5). 
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Substituting these definitions into (16) gives 


— V • cr(xf,(6», s, t), t) = L(u, v) + pA^w(6», s, t) + f(s, t), 
p 'at 


G iP, s G 7^. (17) 

Finally, using (10) for the force f yields the continuous compatibility interface condition 

0 G s&n, (18) 


(crn)(^,s,t) dO —V ■ (T(xfc(^, 5, t), i) = L(u, v) -\- pA^^{9,s,t), 
-•> p ^ ' at 


which can be interpreted as a mixed Robin-type condition along the surface of the beam coupling the 
behaviors of the fluid and structure. 


3.2. Partitioned Schemes Using the AMP Condition 

Following the approach used in [5], we employ an AMP interface condition based on the compatibility 
condition (18) as a stable condition for partitioned schemes coupling incompressible fluids and thin structural 
beams. Equation (18) appears to be a condition that fully couples the fluid and solid values on the interface. 
However, a critical observation made in [5] is that (18) can be used as a stable condition for partitioned 
schemes by using predicted values for the solid variables and when evaluating the right-hand- 

side of (18). These predicted values for the solid are obtained in a first step of the partitioned scheme and 
the interface condition then becomes a mixed Robin-type condition for the fluid in a subsequent step of the 
scheme. This leads to the primary AMP interface condition: 

Condition 1. The AMP interface condition for the fluid is 

[ {(Tn){9, s,t) dO +—V ■ (T{xb{0, sfl),t) = {9, s,t), 9 gP, sgP, (19) 

Jp P at 

where v(p) and are predicted solid variables, and . 

The essential philosophy in this approach is to embed the local behavior of solutions of the FSI problem 
to define the interface dynamics. Therefore, in a partitioned scheme the solid would be advanced first to a 
new time step, providing the predicted values used in the interface conditions (19) for the subsequent fluid 
solve. This process can be viewed as computing the solution of a generalized fluid-solid Riemann problem 
at the interface, thus extending the approach used in the AMP scheme for elastic solids coupled to inviscid 
compressible fluids discussed in [6, 7, 24]. In these previous studies, the FSI problems were purely hyperbolic 
and the fluid-solid Riemann problem involved local, piecewise constant states of the fluid and solid on either 
side of the interface. The solution was then computed using the method of characteristics, giving interface 
values that are impedance-weighted averages of the fluid and solid states. In the present case, involving 
an incompressible fluid, the initial states for the solid domain are and while the solution of 

the generalized fluid-solid Riemann problem requires the solution of an elliptic problem in the entire fluid 
domain with (19) providing the appropriate condition along the interface. 

If a fractional-step scheme is used to integrate the fluid equations, with the velocity and pressure being 
updated in separate stages, then suitable boundary conditions can be derived by decomposing (19) into 
normal and tangential components. The normal component gives a boundary condition for the pressure, 
while the tangential components, along with the continuity equation (2), provide boundary conditions for 
the velocity. 

Condition 2. The AMP interface condition for the pressure equation (4) is given by 

n”^ [ {pn){9,s,t)d9 +—^ = 

Jp P on 

for 9 G V and s G TZ, where 

{pn){9, s, t) = p{iib{9, s, t),t) n{xb{9, s, t), t), ('rn)(0, s, t) = T{xb{9, s, t), t) n{xb{9, s, t),t). 


at 


(p) 


ppA 


Av - 


(Tn)(0, s, t) d9 


(20) 
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Condition 3. The AMP interface conditions for the velocity when integrating the momentum equation (1) 
are given by 


(Tn)(0, s,t) dO + 

’ P 


= tT„ — Vp + + [ {pn){9,s,t) d6 

p at Jr, 


( 21 ) 


for 9 G V and s G TZ, where = tmi9,s,t), m = 2, denote mutually orthonormal tangent vectors. 

We note that if the beam motion only depends on one degree of freedom (e.g. as with the Euler-Bernoulli 
beam model described in Section 5) then only (20) is used and (21) should be replaced with a standard 
condition on the tangential velocity. For example, (21) can be replaced by setting the tangential components 
of the fluid velocity equal to the tangential components of the velocity on the beam surface, i.e. a no-slip 
type condition on the tangential velocity, 

t^v(xh(6',s,t),t) = t^v(s,t)+t^w(6»,s,<), m = l,2, (22) 


for 9 G V and s G TZ. Another choice for the boundary condition could be a slip wall condition which sets 
the tangential components of the traction to zero, 

tm('rn)(6», s, t) = 0, m = 1,2, (23) 


for 9 gV and s G TZ. 


3.3. Velocity projection 

The velocity at the interface in the AMP scheme is defined to be a weighted average of the predicted 
velocities from the fluid and solid. This averaging can be motivated by noting that for heavy beams the 
motion of the interface is primarily determined from the beam and thus the interface velocity should be 
nearly equal to that of the structure. For very light beams, on the other hand, the fluid interface approaches 
a free surface where the motion of the interface is primarily determined by the fluid itself. The weighted 
average given in [5] correctly accounts for these two limits, and this average is extended here to account for 
beams of finite thickness. 


Velocity Projection Algorithm. 

1. Compute the projected surface velocity v^{9,s,t) on the surface of the beam as a weighted average of 


the predicted surface velocities of the fluid vl^\9, s,t) and the beam vl^\9,s,t), 

v^( 6 »,s,t) = 7 v[^^( 6 »,s,t) + (1 - 7 )v[^^( 6 »,s,t), 9gP, sg1Z, 

1 

7 = 


where 


(24) 


(25) 


1 + ipA)/{pAf) 

Here, Af = Af{s,t) is the area of a characteristic annular region in the fluid surrounding the beam at 
a position Xo(s,t) on the reference curve. It has been found that algorithm is insensitive to the choice 
of this characteristic area [5]. 

2. Transfer the interface velocities (minus the surface correction w) onto the beam reference curve 
using an average of the values on the cross-section perimeter curve V 


v.(.,i) = T 


(v\9, s,t) — w{9, s,t)^ d9, 


S 7 ^, 


(26) 


where \'P\ denotes the length of the perimeter curve V. We note that if the fluid density is not constant, 
or if the beam separates two fluid regions in two dimensions, then the average in (26) could be replaced 
by a density-weighted average. Redefine the reference-curve velocity "visfi) to equal Vo(s,t). 








3. Set the fluid velocity on the surface of the beam to he 

v(xfc(0, s, t), i) = v(s, t) + w(0, s, i), 9 s G 1Z, (27) 

following the kinematic matching conditions in (8) and (9) so that the beam reference-curve velocity 
provides a consistent definition of the fluid surface velocity. 

From the perspective of numerical stability, this velocity projection appears to be needed only when 
the beam is quite light. Otherwise the fluid velocity on the interface could be taken as the velocity of the 
beam. This observation is supported by the linear stability analysis in [5] as well as numerical experiments. 
However, in practice we find that the velocity projection is never harmful and we therefore advocate its use. 
In addition, for very light beams the interface acts as a free surface with no surface tension (assuming L is 
also very small). In this case, the fluid interface velocity is smoothed using a few steps of a fourth-order filter 
to remove oscillations that can be excited on the interface, see Section 7.3 for more details. Finally note 
that for cases where the beam supports motion only in the normal direction to the beam, only the normal 
component of the fluid velocity is projected, i.e. 

(n^v)^( 6 »,s,t) = 7 (n'^v[^^)( 6 »,s,f)-h (1 - 7 )(n^v[^^)( 6 »,s,t), OgV, sgP, (28) 

where 7 is given by (25) as before. 


4. The AMP FSI time-stepping algorithm 


In this section we outline the FSI time-stepping algorithm, including deforming geometry, for an incom¬ 
pressible fluid governed by the Navier-Stokes equations and a thin structure determined by the general beam 
model in (5). The intent is to provide a concrete description of the AMP algorithm with sufficient details 
to indicate the primary extensions required from the original scheme described in [5]. The fluid equations 
are solved in velocity-pressure form using a second-order accurate fractional-step algorithm [5, 23, 25]. The 
beam equations are advanced with a second-order accurate predictor-corrector scheme. In practice we have 
also used a predictor-corrector version of the Newmark-beta scheme for the beam equations. 

Given the discrete solution at the current time t", and one previous time level the goal of the FSI 
algorithm is to determine the solution at time The algorithm is written for a fixed time-step. At, so 

that t” = nAt, but is easily extended to a variable At. A predictor-corrector scheme that implements the 
FSI algorithm is as follows: 

Begin predictor. 

Stage I - structure: Predicted values for the displacement and velocity of the shell are obtained using 
a second-order accurate leap-frog scheme given by 


-Ap) 




2At 


j G 


' 2At' = j ^ 

where the external force at t = is given in terms of the fluid stress on the surface of the beam by a 
discrete approximation to ( 10 ), 


fj"= {rTn)ZA0. 

k G Vh,] 

Here j G TZh denotes the indices of the discrete points along the beam reference curve and i G F/i denotes 
the corresponding indices of the discrete fluid points on the beam surface. There is a cross section of the 
beam identified for each point j along the reference curve, and the points k G Vhj denote the discrete 
representation of the curve on the perimeter of the cross section. 
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Stage II - grid generation: Predicted values for the deformed interface positions are obtained based 
on the predicted reference-curve deformation using discrete versions of the formulas in ( 6 ) and (7). The 
predicted fluid grid is then defined in terms of the interface positions using an appropriate mapping. 
In the context of overlapping grids (see Section 7.1) an interface-fitted fluid grid is generated with a 
hyperbolic marching algorithm and then a new overlapping grid is generated. 

Stage III - fluid velocity: Predicted values for the fluid velocity are obtained using a second-order 
accurate Adams-Bashforth scheme 


vp) _ vP 
^ At 


Q 1 

y-pin _ _Tpn—1 

2 ' 2 ' 


i G 


where 


FP = -((uP - wP) • V)uP - V„pP + /iA^^vP. 


On moving grids the equations are solved in a moving coordinate system (see Section 7.1) and wP is the 
grid velocity. From (21), the boundary conditions for on P/i are 


^m,i 

^ r(f)nkA0-f 

- t"^ . 





L ^ kGT^.j J 


V„ • = 0, 


where i G P^,, j G TZh, m = 1 , 2 , and 


(P) 

T; ^ M 






- (P) 
■■'V. ^ 


(29) 


(30) 


The term in (30) is a predicted value for dw/dt, which for an Euler-Bernoulli beam can be computed 
from the predicted beam position, velocity and acceleration. At this stage of the time-stepping algorithm, 
a predicted value for the pressure, p\^\ on the right-hand side of (29) is not yet available. However, a 
sufficiently accurate value for the purposes of the boundary conditions can be obtained using a third- 
order extrapolation in time, = 3pP — 3pP“^ +Pi~^- This extrapolated pressure is replaced in the 
next stage of the algorithm by solving a discrete Poisson problem. Appropriate boundary conditions for 
V; should also be applied on other boundaries. 

Stage IV - fluid pressure: Predicted values for the pressure are determined by solving the pressure 
equation 


= -pVvp^ : , ieflh, (31) 

subject to the boundary conditions. 


T 

nf 

Pk ^nk A6» -H — 

T 

= ni 



^ k G Vh,i ^ ^ 


^ ^ k G Vh,i ^ 


with j G TZh and appropriate conditions for p: applied on the other boundaries. Note that in practice 
it is useful to add a divergence damping term to the right-hand side of the the pressure equation (31) 
following [23, 25]. 

End predictor. 

Begin corrector. 
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Stage V - structure: Corrected values for the displacement and velocity of the shell are obtained using 
a second-order accurate Adams-Moulton scheme, 


pA 


u!*+i - Uj” 
At 

v!"+i - 

* ' At ' 


-"+5 

= ^j ’ 


= L, 


/_n+i 

h ’^j ) 




j G 

j G T^ht 


where 


-”+5 

"j 


5 (“]'’+<) 


-"+5 

''j 


(v“ + v“) , 


3^-1-i 

fj = 


1 (fip) 

2 V J 




Stage VI - grid generation: The fluid grid could be corrected and regenerated at this stage of 
the algorithm. In practice, however, this stage is generally skipped to avoid the cost of an extra grid 
generation step. Note that the predicted grid is already second-order accurate and there appears to be 
no noticeable adverse effect on the stability of the algorithm by skipping this stage. 

Stage VII - flnid velocity: Corrected values for the fluid velocity are obtained using a second-order 
accurate Adams-Moulton scheme, 


vr+^ - vr 


At 


2 ^ * 


■FD, 


i G flh. 


The boundary conditions have the same form as those used in Stage III with uj^^, and replaced 
by u"'*'^, and Wj"^ 

Stage VIII - fluid pressure: Corrected values for the pressure are determined by solving the discrete 
equations in Stage IV with the predicted values replaced by values at 

Stage IX - project interface velocity: Define an interface velocity as a weighted average of the fluid 
and solid interface velocities, using 

= yvl'+l + (1 - 7) (v;+l + w;+l) , iGT;,, 


where 7 is defined in (25). Redefine the beam reference-curve velocity as an average of the surface values 
(adjusted for 


v”+i 


I^Aj 


E 




A0, 


j G TZh, 


k G Vh,j 

and then redefine the fluid surface velocity to be consistent with the new beam velocity, 

iGT/i, iGTZh- 


v”+^ = Vj”+^ 


■w;+\ 


End corrector. 

We emphasize that the AMP algorithm is stable with no corrector step, although if the predictor step is 
used alone, then Stage IX should be performed following the predictor to project the interface velocity. We 
typically use the corrector step, since for the fluid in isolation the scheme has a larger stability region than 
the predictor step alone, and the stability region includes the imaginary axis so that the scheme can be used 
for inviscid problems {p = 0). 
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5. An Euler-Bernoulli beam model 


For the purposes of this paper, we consider a two-dimensional beam governed by a generalized Euler- 
Bernoulli (EB) beam model. The model considers the behavior of the displacement, ri{s,t), of the reference 
curve of the beam in the direction normal to its initially flat position, where s measures distance along 
the reference curve. The model assumes small deviations from its initial position (i.e. small slopes) so that 
displacement in the direction tangent to the reference curve is neglected. The generalized EB model includes 
damping terms and takes the form 


ph 


d'^T] 

W 


Kop + 





as2 ) 





as2 





s e Tl, (32) 


where p is the density of the beam, h = h{s) is its thickness, Kq is the linear stiffness coefficient, T is the 
tension coefficient, E is Young’s modulus, and I is the area-moment of inertial of the beam. The term with 
coefficient Ki is a linear damping term, while the term with coefficient Ti damps higher spatial frequencies 
in ??(s,t) faster than lower frequencies and thus serves to smooth the beam. We assume that all of the 
coefficients in (32) are constants, except for the thickness of the beam which may vary with s such as near 
the rounded end of a cantilevered beam (see Figure 5, for example). 


x+(s,i) 

xo(s,t) 

x_(s,i) 


Figure 3: Geometry of an Euler-Bernoulli beam. The beam reference-line is xq = xo(s,f) and the beam normal is n. 
The beam has thickness h and meets the fluid on the curves x_ and X 4 -. Outward normals to the fluid domain are 
n_ and n+. 



The geometry of an FSI problem involving a section of an EB beam is illustrated in Eigure 3. For this 
configuration, the displacement of the reference curve is given by u(s,t) = (0,r]{s,t)) so that its position is 
xo(s,t) = (s, ry(s,t)). The unit normal to the reference curve, n(s,t), can be computed in terms of the slope 
drj/ds, and this determines the position of the upper and lower surfaces of the beam as 


x±(s,t) = Xo(s,t) -bz±(s,t), z±{s,t) = ±^n{s,t), sell, (33) 

which is a specific case of the formula given previously in (12) with the b subscript dropped for notational 
convenience. From (33), the velocity of the upper and lower surfaces of the beam are given by 


d h d 

v±(s,t) = —x±(s,t) = v(s,t)-b w±(s,t), w±(s,t) = ±-—n(s,t). (34) 

Finally, a modeling choice is required to determine the force /(s,i) on the reference curve of the beam, and 
we use 




(n’^crn)^ -b (n'^crn)^ 


(35) 
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where n is the unit fluid normal on the beam surface so that (n^crn)± is the component of the fluid traction 
in the direction normal to the upper and lower surfaces of the beam. This modeling approximation is valid 
for small beam deflections, and it is consistent with the use of an EB beam model and (a reduced form of) 
the AMP pressure condition in (20) where the right-hand side of (35) becomes the left-hand side of (20). 

For a cantilevered beam, the free end of the beam is usually rounded as shown in Figure 5. The flow 
can be complex near a sharp corner and by rounding the corner a high resolution grid can be generated that 
resolves the flow in this region. In addition, the force on the beam surface is computed using the beam-surface 
normal n and the approximation n^crn. Near the rounded end of a beam this is a poor approximation since 
the fluid normal may be orthogonal to the beam reference curve. Consequently, in practice, the fluid forces 
on the beam tip are ignored in computing force on the beam (in particular, any points on the surface of 
the beam that extend beyond the length of the beam reference line are ignored). This should be considered 
a modeling assumption. Similarly, contributions from the end points are also ignored when projecting the 
interface velocity onto the beam. 

6. Analysis of a model problem for a two-dimensional beam 



Figure 4: The geometry of the linearized model for a beam with fluid on two sides. 

In this section, the stability of the AMP algorithm is investigated for a linearized model problem con¬ 
sisting of an EB beam coupled to inviscid, incompressible fluids with different constant fluid densities in the 
regions above and below the beam as shown in Figure 4. For comparative purposes, we also consider the 
stability of a traditional partitioned (TP) scheme for the same model problem. We begin by obtaining the 
exact solution of the model problem to explicitly identify the contribution of the added mass of the fluid 
for two different choices of the boundary conditions on the fluid. The stability of the discrete AMP and TP 
schemes is then investigated using a mode analysis that extends the theory for the case of a beam with fluid 
on one side only presented in [5]. For both the AMP and TP schemes, the fluid equations are advanced using 
a backward-Euler time integrator as opposed to the predictor-corrector integrator described in Section 3. 
While this simplifies the analysis somewhat, the essential results are unaffected. It is found, in agreement 
with [5], that the AMP scheme is stable for arbitrarily light structures when the usual time-step restriction 
is satisfied. The TP scheme without sub-time step iteration, on the other hand, becomes unconditionally 
unstable for light structures. Note that the TP scheme has been analyzed previously by Causin, Gerbeau 
and Nobile [1], among others, and the results presented here extend those works. 
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In the following analysis the two fluids are coupled to an EB beam with constant thickness h and without 
damping terms. The equations describing the evolution of the fluid-structure system are then linearized about 
flat undeflected beam surfaces at j/ = —h/2 and h/2 as shown in Figure 4. Let v+(a::, y, t) and p'^{x, y, t), 
denote the fluid density, velocity and pressure, respectively, in the upper domain, 11 + = [ 0 , 1] x [h/2, H^], and 
p~, w~{x,y,t) and p~{x,y,t), the corresponding quantities in the lower domain = [0,/] x [H~, —h/2], 
where I is the length of the beam. The governing equations and boundary conditions for this model are 


Fluid: 


Beam: 

Interface: 


V • v± = 0, 

= 0 or = 0 , 


X e n±, 

X e n±, 

X G (0,1), y = 


= -^oV + f^-El^-[p]^, xG{0,I), 


dx"^ 


dr] 


"^ 2 =^^ a;e( 0 ,Z), y = ±'- 


dx^ 
h 
2 ’ 


(36) 


where 


[p]p V 2 ,t)-P {x,-h/2,t), 

denotes the jump of the pressure across the beam, and where we take periodic boundary conditions in the 
^-direction. Initial conditions for v^, 77 and dp/dt are needed to complete the specification of the problem, 
but these are not required in the subsequent analysis. 


6.1. Exact solution of the continuous model problem 

It is instructive to first examine the behavior of the exact solution of the FSI problem in (36) . Since this 
model problem is a linear constant-coefficient system of equations with periodic boundary conditions in the 
a;-direction, the solution can be expressed in terms of a Fourier series in x, e.g. 

00 

^^{x,y,t)= ^ v=^(fc,y,t)e*'""^, kx = 2Trk/I, 

k— — (x> 


and similar expressions for p^ and 77 ^. Here, kx identifies the wave number of the Fourier mode. Transforming 
the equations for the fluid to Fourier space gives 


dt 




■7 -± 

ikxvf + 


dp^ 

dy 


= 0 


dy 


= 0 


y G {H , —h/2) for the lower fluid, 
y G (h/2, iL+) for the upper fluid. 


(37) 


The transformed equation for the beam is 


ph 


d'^f) 

W 


-if]- 


(38) 


where L = Kq + Tk^ + Elk'/., and for notational brevity we have suppressed the functional dependence on 
the wave number. The transformed equations for the fluid in (37) can be manipulated to derive equations 
for the pressures of the form 


d'^p'^ 

dy"^ 


S,2+± 


= 0 . 


(39) 
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It is assumed that kx ^ 0 in the present analysis while the special case when kx — 0 is considered later. 
For the case of the rigid-wall boundary conditions, = 0 at y = we use the y-momentum equations 
in (37) to derive the Neumann conditions for the pressures, 




(40) 


General solutions of the pressure equations in (39) satisfying the homogeneous Neumann boundary conditions 
in (40) are 

cosh {kxiy - H^)) , (41) 

where = a^{t) are a time-dependent coefficients. Similarly, the kinematic conditions at the interface, 

vf{±h/2,t) = 

together with the y-momentum equations in (37) give the interface conditions 


(42) 


The general solutions in (41) and the interface conditions in (42) imply 

a± = csch (kxV^) , 


(43) 


where = zLH^ — h/2 measure the vertical extent of the upper and lower fluid domains. Using (41) 
and (43), the applied force appearing in the beam equation (38) takes the form 


Titr ntr-l- 


where 


Mt = 


“ dt^ 
coth{kx'D^) 


kxV± 


dt^ 


(rigid-wall BCs), 


(44) 


(45) 


are the added masses from the upper and lower fluids for the case of rigid-wall boundary conditions at 
y = H^. Note that represents the total mass of fluid in the upper and lower domains. For the 

alternate case of the pressure boundary conditions, = 0 at y = , a similar analysis leads to the added 

masses given by 




tanh(fca;I?^) 


(pressure BCs). 


kxV^ 

Finally, using the applied force (44) in the equation of motion for the beam (38) gives 


{ph 


+ M- + M+ 


\ _ 
) df^ 


Lfj, 


(46) 


(47) 


which explicitly identifies the contribution of the added mass of the fluid in the beam equation. A solution 
to (47) satisfying initial conditions for t) and dp/dt can be found easily to complete the solution for r). Given 
r), the pressure is then fully determined from (41) and (43), which in turn determines and from (37) 
together with the initial conditions and boundary conditions. This completes the solution of the model 
problem in Fourier space. 

We now return to the special case when kx = 0. For rigid-wall boundary conditions, the full solution for 
fcr = 0 is found to be 


= ^f(y,o). 


V2 = 0 . 


?7 = y(0), p±=p±(t), 


(48) 


where pj (t) are spatially uniform fluid pressures whose difference Apo = Pq" — Pq is a constant related to 
the beam displacement by 

KoviO) = -Apo. 
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Note that the kx = 0 mode of the beam is stationary, consistent with the fact that the volumes of the 
lower and upper domains must remain constant for incompressible fluids. The added mass in this case is 
interpreted to be infinite. This is fully consistent with the observation that the added mass in (45) tends to 
infinity as kx'D^ tends to zero. For the case of pressure boundary conditions with kx = 0, solutions for the 
fluid pressures become 



so that the added masses in (44) are 

This is in agreement with the limit kx'D^ 0 for the formula for in (46). 

Returning to the equation for fj(t,kx) in (47), we note that added mass associated with the fluids on 
either side of the beam tends to reduce the acceleration of the beam. Since —)■ 0 as kxT)^ —>■ oo for 
either (45) or (46), the effect of the added mass is small for high-frequency components of the solution. 
For the case of rigid-wall boundary conditions, —)■ oo as \kx'D\ becomes small so that the effect of the 

added mass is large in this limit. This suggests that the rigid-wall case could be difficult numerically if the 
algorithm is sensitive to the effects of added mass. For the case of the pressure boundary condition, 
becomes constant when kxTP = 0 and so this case would likely be less difficult numerically than the case of 
rigid walls. 

6.2. Stability analysis of the AMP and TP schemes 

The stability of a partitioned scheme for an FSI problem depends crucially on the treatment of the 
interface and the resulting boundary conditions applied to the fluid and structure on either side. AMP 
interface conditions have been developed and used in the AMP FSI time-stepping algorithm, and the purpose 
of this section is to analyze the stability of a form of this algorithm applied to the model problem in (36). A 
starting point for this analysis is the system of equations in (37) and (38) already written in Fourier space. 
For this system, we consider a discretization in time, but let the dependence on y remain continuous. This 
approach reveals the stability of the partitioned algorithm, while keeping the algebraic details as simple as 
possible. 

We begin with an analysis of the AMP algorithm, and then compare the results with a traditional 
partitioned (TP) in the end. For both cases, let 

v-±(2/)«v±(y,r), p"±(y)«p±(y,r), ^ ^ fin, 

where t" = nAt. Furthermore, let D+t and D-t denote the forward and backward divided different operators 
in time (e.g. Dj^tfp = ( 77 ”+^ — fp)/ At and = ( 7 )” — Tp~p/ At). 

The following AMP scheme is an abbreviated version of the full scheme described in Section 4: 
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AMP scheme: 

Stage I: Advance the beam displacement using a leap-frog scheme, 

phD+tD.tT = -Lr -[P"]r- ( 49 ) 


Stage II: Advance the fluid velocity and pressure using a backward-Euler scheme applied to the 
equations in pressure-velocity form, 


Qpin+l)± 
p^D+tV2 + TT 






0 


0 


y G {H ,—hf2) for the lower fluid, 
y G {h/2,H~^) for the upper fluid. 


(50) 


with the AMP Robin interface conditions derived from (20), 


and boundary conditions. 


dp 

dy 


phdp^^+^^^ 

dy 

= Lr+\ y = ±^, 

(51) 

= 0, 

y = H^. 

(52) 


A key ingredient of the AMP scheme is the Robin condition in (51) for the fluid pressure. This condition, 
along with (49), implies 

/ j/ = ±-, (53) 

dy 2 

which is a discrete version of the interface condition in (42). We observe that the acceleration term on the 
right-hand side of (53) is evaluated at the same time level as the pressure gradient term on the left-hand 
side, and this plays an important role in the stability of the AMP scheme as is noted below. After solving 
the pressure equation for the fluid subject to the boundary condition in (52) and interface condition in (53), 
the applied force is found to be 

- ] r = (54) 

where the added masses, are given in (45). Substituting the expression for the applied force (54), with 
n -I- 1 replaced by n, into the leap-frog scheme (49) gives the difference equation for ?)”, 

[ph + M- + M+) D+tD.tr = -Lr, (55) 

which is analogous to the evolution equation for fi{t) in (47). A solution of the form ?)" = A^rj^ is sought 
which leads to the quadratic equation 


- 2 1 - 


AfL 


2{ph-\-Ma 


A-kl = 0, 


(56) 


for the amplitude A. Stability requires |A| < 1 and that the roots be simple, if they have modulus one, 
which leads to the following theorem: 
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Theorem 1. The AMP scheme applied to the model problem in (36) is stable if and only if 


At <2 



( 57 ) 


where is given by (45) for case of the rigid-wall BCs and by (46) for the case of pressure BCs. Moreover, 
if (57) holds, then the roots of (56) satisfy |A| = 1 and the AMP scheme for the model problem is non- 
dissipative. 

We note that if a stable time step is chosen for the leap-frog scheme in (49) with no external forcing 

due to the fluid, namely At < 2\Jph/L, then this value of At would be stable for the AMP scheme as well 

since Mf; > 0. We also note that a velocity projection has been omitted in the prototype AMP scheme 
considered in the stability analysis, and thus the fluid and structure velocities at the interface would only 
match to first-order accuracy. An additional projection step to match the velocities on the interface exactly 
could be included, but this is not essential for the present analysis. 

Having considered the stability of the AMP scheme, we now turn our attention to a TP scheme defined 
as follows: 


Traditional partitioned (TP) scheme: 

Stage I: Advance the beam displacement using (49). 

Stage II: Advance the fluid velocity and pressure using (50) with interface condition, 


dy 


^ = -p^D+tD-tT, 


h 

y = ±^. 


(58) 


and the boundary condition in (52). 


The only difference between the AMP scheme and the TP scheme is the coupling at the interface. For 
the TP and AMP schemes, the stress from the fluid at time level t” provides a forcing to the leap-frog scheme 
in (49) to advance the displacement of the beam to time level The velocity and pressure in the fluid 

are then advanced to time level using the discrete equations in (50), the boundary condition in (52) 
and a coupling condition at the interface. The AMP scheme uses the coupling condition in (51), whereas 
the TP scheme uses (58). The latter is derived from the time derivative of the kinematic interface condition 
matching the vertical velocity of the beam to that of the fluid on interface. We note, however, that a similar 
interface condition was derived from the leap-frog scheme in (49) and the AMP interface coupling condition 
in (51), which appears in (53). The key difference is that the time level for the vertical acceleration of the 
beam in (58) lags that of fluid in contrast to (53), and this is the essential issue leading to the instability of 
the TP scheme for light structures as the subsequent analysis shows. 

Following the steps used in the analysis of the AMP scheme, the pressure equation in (50) is solved with 
interface conditions in (58) and boundary conditions in (52) to give the applied force 

- ] r = -M-D+tD.tr - M+D+tD_tr- (59) 

Using the beam equation in (49) together with (59) to eliminate the applied force yields a difference equation 
for ?)” of the form 

phD+tD_,r + (m-+ M+) D+tD_tr~^ = -ir- (eo) 

Again we observe a time lag in the added-mass terms in (60) as compared to the previous evolution equation 
for ff^ in (55) for the AMP scheme. This time-lag is well known and is the primary source of the added-mass 
instability for the TP scheme. Using r)" = A'^ff in (60) leads to the cubic equation for A, 

phA{A - if + At^LA^ -f {M- + MffA - if = 0. (61) 
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The TP scheme is weakly stable if all roots of (61) satisfy |^| < 1, and with no repeated roots of modulus 
one. Using the theory of von Neumann polynomials [26, 27] leads to the following theorem. 

Theorem 2. For ph > 0 and L > 0, the traditional partitioned scheme applied to the model problem in (36) 
is weakly stable if and only if 

M- + M+ < ph (62) 


and 


At < 2\ 


Iph-MF -Mf 


(63) 


where is given by (45) for case of the rigid-wall BCs and by (46) for the case of pressure BCs. 

The condition in (62) implies that the TP scheme is unconditionally unstable if the contribution from 
the added mass is too large, i.e. if M~ + M+ > ph, regardless of the choice for the time-step. At. In 
general, the results of the analysis suggest that TP-type schemes would suffer from stability issues when 
M^/ph is large. We note from (45) or (46) that Mf; /ph is large if the mass of the fluid domain is large 
compared to that of the structure, i.e. if p^V^/ph is large. However, the ratio is also large for the case of 
rigid-wall boundary conditions if \k,j;'D^\ is small in view of the bracketed term in (45). This latter condition 
can be re-interpreted as V^/l is small, i.e. the fluid domain is long and slender, since kx = 2'Kk/l. These 
observations are in agreement with the results of Causin et al. [1], who studied a simple beam with fluid on 
one side. The observations are also in agreement with the numerical results presented here and in [5]. If the 
condition M~ M+ < ph is satisfied, then (63) indicates how the time-step must be reduced as the total 
added-mass increases. 


6.3. Pressure regularization for rigid-wall boundary conditions 

Note that for both the TP and AMP schemes the limit kx = 0 with rigid-wall boundary conditions 
must be treated with some care since the continuous fluid pressure is only determined up to an arbitrary 
additive constant. In the AMP scheme p'^ and p~ are coupled so that only the difference in pressure is 
arbitrary, whereas both pressures are arbitrary for the TP scheme. In a numerical implementation these 
singularities need to be removed since, due to round-off or truncation errors, the singular systems may have 
no solution. The specific choice of regularization can affect the overall stability of the approach. For example, 
one common regularization would follow the approach in [25], which for kx = 0 would solve the augmented 
pressure equations 

g2^(n+l)± 

dy2 + “ 


with the constraints 


P 


(n-Hl)- 


^H+ 


dy = 


= 0. 


IH- 


Here are two additional unknowns which serve to regularize the solution and ensure that the aug¬ 

mented system has a unique solution. The effect of regularization can be significant for the kx = 0 mode. 
Using this regularization the traditional partitioned scheme is found to be stable for the kx = 0 mode under 
the mollified time-step restriction 


At <2 



This is significant difference from (63) which required At —>■ 0 when kx —>■ 0. This less severe time-step 
restriction is a result of the pressure regularization and applies only to the kx = 0 mode. All other small 
kx modes in the system must satisfy (63). An analogous regularization for the AMP scheme involves the 
pressure equations 


^2^(n-|-l)± 


+ Q-("+l) = 0, 
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and the single constraint 


P 


("+!)- 


^H+ 


dy + 


P 


(n 


+i)+d?/ = 0, 


fH- 


which determines the one unknown In this case it can be shown that the AMP scheme is stable with 

no time-step restriction for = 0, and the solution is found to be 


.(n-|-l)± 


r-.Oi 


d"+l)± 


= 0 , 


fin+1 ^ 


P 


(-+i)± = 




V++V- 


Note that, by adding one pressure regularization, the AMP scheme is able to determine the solution uniquely, 
and the solution is compatible with the exact solution for kx = 0 given in the continuous case by (48). 


7. Numerical approach using deforming composite grids (DCG) 

Our numerical approach for the solution of the equations governing an FSI initial-boundary-value problem 
is based on the use of deforming composite grids (DCG). This FSI-DCG approach was first described in [6] 
for the case of an inviscid compressible flow coupled to a linearly elastic solid, and later in [7] for the case of 
compressible flow coupled to nonlinear hyperelastic solids. Here, we extend the approach to FSI problems 
involving an incompressible flow coupled to deforming beams. 

7.1. Deforming composite grids and the fluid domain solver 

Deforming composite grids (DCGs) are used to discretize the evolving fluid domains in physical space. 
An overlapping grid, Q, consists of a set of structured component grids, {G^}, g = 1,... ,JV, that cover each 
fluid domain, D/;(t), and overlap where the component grids meet. Typically, boundary-fitted curvilinear 
grids are used near the boundaries while one or more background Cartesian grids are used to handle the bulk 
of the fluid domain. Each component grid is a logically rectangular, curvilinear grid in space dimensions, 
and is defined by a smooth mapping from parameter space r (the unit square or cube) to physical space x, 

x = g(r,t), rG[0,l]”G xGK"A 

Typically, the background grids are static, while the boundary-fitted grids evolve in time to match the motion 
of the boundary. 



Figure 5: Composite grids for a beam in a channel at times t = 0.0, and 1.0. The green fluid grid deforms over time 
to match the evolving beam (shown in white) and overlaps with the background Cartesian fluid grid (shown in blue). 
The reference curve of the beam is shown in red and the thickness of the beam describes the position of the surface 
of the beam. 

In the FSI-DCG approach, component grids next to a fluid-structure interface deform over time to match 
the beam motion. This is illustrated in Figure 5 for the case of a beam in a fluid channel. A green fluid 
grid is fitted to the boundary of the beam, which is defined in terms of displacement of the reference curve 
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of the beam (shown in red) and the thickness of the beam (c.f. (33)). The beam thickness is a function 
of the position s along the reference curve, and remains fixed as the beam deforms. After each time step, 
the points on the fluid interface corresponding to the surface of the beam are recomputed according to the 
current displacement of the beam and the normal to the reference curve using (33), and a hyperbolic grid 
generator [28] is used to regenerate the local interface grids. After all deforming component grids have 
been regenerated, the Ogen grid generator [29] is called to regenerate the overlapping grid connectivity (e.g., 
cut holes, determine interpolation points) between the evolving interface-fitted fluid grids and the fixed 
background fluid grids, typically Cartesian grids, shown in blue in Figure 5. 

The incompressible Navier-Stokes equations are solved in velocity-pressure form on overlapping grids 
using a second-order accurate fractional-step algorithm [5, 23, 25] similar to that described in Section 4. 
The treatment of moving overlapping grids, including the assignment of exposed points, follows the approach 
discussed in [30] for compressible flow. The governing equations for the fluid on a given component grid are 
solved in a coordinate frame moving with the grid. Consider the momentum equations in (1) for the fluid in 
an Eulerian frame. 


dvi dvi 1 daji 

dt dxj p dxj ’ 


i 1 ,..., Tifi, 


(64) 


with summation convention. Under a general moving coordinate transformation, x = g(r, t), the equations 
in (64) are transformed using the chain rule to become 


dvj 

iH 


+ 


Wj) 


dvk dvi 
dxj dxk 


1 dvk daji 
p dxj dxk ’ 


i 1 , . . . , TLfl, 


(65) 


where the component of velocity Vi and stress Uji are now considered to be functions of r and t, and where 
Wj = dgj/dt are the components of the grid velocity, w. The pressure equation in (4) is transformed to (r,t) 
in a similar way. 


7.2. Structural domain solvers and the AMP condition 

Two approaches are implemented to advance the governing equations for the beam. The first is a 
second-order accurate predictor-corrector scheme described previously in Section 4, which also uses standard 
second-order accurate finite differences to approximate the spatial derivatives in (32) assuming an Euler- 
Bernoulli beam. The basic AMP condition in (19) for this beam model becomes 

(crn)+ 4- (crn)_ 4-— V • (T(x±(s,t),t) = L(u(P\ -I- p/i^w2’^(s,t), sgH, (66) 

p ^ ' at 

where the predicted quantities on the right-hand side of (66) involving the displacement and velocity of 
the beam can be computed in terms of the discrete values for j], drj/dt, and d^rj/dP, and their spatial 
derivatives, using standard finite differences. The normal and tangential components of (66) then become 
boundary conditions for the fractional-step fluid solver. 

The second approach implemented to advance the equations in (32) for an Euler-Bernoulli beam uses a 
standard finite-element approximation in space and the Newmark-beta time-stepping scheme (with param¬ 
eters chosen for second-order accuracy in time). The spatial representation of the beam displacement in the 
finite-element approximation is taken as 


Vh{s,t) = }, (67) 

i=i 

where the degrees of freedom are the nodal displacement r]j{t) and slope pi it), and where (('j(s) and '4’jis) 
are cubic Hermite polynomials with support on the interval [sj_i, s^+i] and satisfying 

Aisj) = 0 , 
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= 0 , 

V'z('Sj) = dij. 




where Sij is the Kronecker delta. Following the usual Galerkin approximation in space, the beam equation 
become 


= -KhVh - ( 68 ) 

where Mfi is a mass matrix, is a stiffness matrix, and Bfi is a damping matrix. The vector contains 
the nodal displacements and slopes, r]j(t) and j = 1,, Nn, and the components of the forcing vector 

ih are given by integrals of /(s, t) in (32) with the Hermite basis functions®. The ODEs in (68) are advanced 
in time using a second-order accurate Newmark-beta scheme. 

Given predicted values for and /dt from the first stage of the AMP algorithm, pointwise values 
of the beam operator, are required to evaluate the right-hand side of the AMP condition in (66). This 
is generally a straightforward computation in the context of finite difference or finite volume approximations 
to the beam operator, whereas the calculation may be more difficult within the finite-element approximation. 
As an example, consider a beam operator of the form 


A — (A/7/5s)ss, 


where the subscripts denote partial derivatives with respect to s. It is sufficient to only consider the highest 
derivative term in (32) for the purposes of the present discussion. Let L have the Hermite finite-element 
representation 

Nn 

Lh{s,t) = ^ }. (69) 

i=i 

where 7 j(t) and 7j(t) are coefficient functions of time, and ipjis) and ^/’j(s) are the cubic polynomial basis 
functions in (67). The 2A^„ unknowns, 7 j(t) and 7j(t), are determined from the 2Nn Galerkin equations 

{Xj, Lh).j^ = iXj,-{Eir]h,ss)ss)Ti, J = 1,2, ...,A^n, (70) 

where Xj = 4'j ot and where (/, g)^^ = /g /(s)g(s) ds. Integration by parts (twice) gives 


{Xj,Eh).^ = - iXj,ss,EI'qh,ss)Ti + { - XjiEIr]h,ss)s + Xj,sEIr]h^ss)'^ 


J = l,2,...,Nr,. 


(71) 


The problem now arises in computing second-order accurate approximations to the boundary terms 
{EIrjh,ss)s and Elrjh.ss appearing in (71). The cubic Hermite representation is not of sufficient accuracy for 
this purpose since, for example, "qh.sss is constant over an element and at best first-order accurate. 

There are various possible solutions to the difficulty of evaluating the right-hand side of (66) in the 
context of a finite-element approximation, and we have considered and tested two different approaches that 
require relatively minor changes. The properties of the two approaches are evaluated in Section 8 where 
some conclusions are drawn. 


7.2.1. AMP-PBA : Predict Beam Acceleration 

In the first approach, which we refer to as AMP-PBA (Predict Beam Acceleration), the vector form 
of (32) is used in (66) to eliminate which gives an AMP condition of the form, 

((Tn)+ + ((Tn)_ -b ^ V • cr(x± (s, t),t) = - f (s, t) + ph^w^^^ (s, t), s^U. (72) 

The predicted beam acceleration jn (72) is readily computed from the acceleration term in the 

weak form (68), while the predicted fluid forcing, f*-P^(s,t) = — (crn)(|f^ — ((Tn)(f\ can be obtained from the 
discrete fluid stress. This gives the following AMP-PBA algorithm: 


“Appropriate adjustments are made to the matrices and right-hand-side to account for boundary conditions. 
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1. Obtain predicted values for = — ((Tn)^|f^ — (crn)^^ using extrapolation in time during the predictor 
stage or the current guess during the corrector stage (c.f. the algorithm in Section 4). 

2. Solve ( 68 ) for /dt^ using the predicted force and evaluate pointwise values for ph{d^\i 
from the Hermite representation for this quantity. 

3. Evaluate the right-hand side to (72) and use this when applying the AMP conditions to the fluid. 

7.2.2. AMP-PBF : Predict Beam internal Force 

In the second approach, referred to as AMP-PBF (Predict Beam internal Force), a mixed finite- 
element/finite-difference approach is used to compute the beam internal force directly from equation (71) 
(extended to include all terms in the full EB model). In particular, using the solution nodal values and 
derivatives on the boundary and two adjacent points, a fourth-order accurate approximation to Elrjss and 
a second-order accurate approximation to {EIrjss)s on the boundary can be derived. This allows the coeffi¬ 
cients 7 j(t) and 7 j (t) to be determined and then pointwise values of can be computed from (69). This 
gives the following AMP-PBF algorithm: 

1. Given predicted values for rj]^\ determine second-order accurate approximations to the boundary 
terms {EIrjh,ss)s and Elrjh.ss in (71) using a finite-difference approximation involving the boundary 
points and two adjacent points. 

2. Compute values for by solving the Galerkin equations in (71) (extended to include all terms in 
the full EB model) for 7 j and 7 '. 

3. Evaluate the right-hand side to ( 66 ) using (69), and use this when applying the AMP conditions to 
the fluid. 

7.3. Filtering the projected interface velocity 

The projected interface velocity (24) is optionally smoothed using a few iterations of a fourth-order 
filter as described in [31]. The beam solution can also be smoothed to prevent high-frequency numerical 
oscillations, for example by choosing the coefficient Ti in (32) to be proportional to As^. However, for 
light beams the interface velocity is determined primarily from the fluid and thus the filter on the interface 
velocity is needed for light beams to prevent numerical oscillations on the interface. Note that as the beam 
becomes very light, with beam parameters ph, T, El, etc. going to zero, the AMP interface condition 
approaches that of a free surface boundary condition for the fluid. Such a free surface condition with 
no surface tension is susceptible to physical and numerical instabilities, and the fourth-order filter acts to 
suppress these instabilities and thus serves as a numerical surface tension. This filter does not affect the 
overall second-order accuracy of the AMP ESI time-stepping scheme. 

7 . 4 . Time step determination 

The global time step. At, for the AMP ESI time-stepping scheme is chosen to be the minimum of the 
stable time steps for the fluid domain solver and beam solver. The At for the fluid solver is chosen to be 
a factor Asf (“SF” for stability factor) times an estimate of the largest stable time step based on a von 
Neumann analysis of the linearized Navier-Stokes system. Typically we choose Asf = 0-9- In this work, 
the FEM beam solver is implicit in time and stable for any time step. However, for accuracy reasons it is 
convenient to choose At based on a condition relating to explicit time-stepping. Thus, the time-step for the 
beam is chosen as a factor Asf times an estimate of the largest time step such that an explicit time-stepping 
procedure for the beam would satisfy a GEL stability condition. Here Asf is chosen in the range 1 to 10. 

8. Numerical Results 

We now present the results for a series of simulations chosen to demonstrate the properties of our 
numerical approach for ESI problems based on the use of the AMP interface conditions. In order to illustrate 
the basic properties of the approach in a simple setup, we begin with two cases that use a single deforming fluid 
grid coupled to an Euler-Bernoulli (EB) beam with negligible thickness and finite mass. This simplification 
eliminates the finite-thickness corrections given in (33) and (34) by setting h = 0. The first test involves a 
time-dependent problem in which an exact solution is known following the method of manufactured solutions, 


23 


while the second test considers a steady state FSI problem with known analytical solution. The solutions for 
both tests are obtained using the beam solver based on finite differences, and they illustrate the accuracy 
and stability of our numerical approach for a relatively simple geometric configuration. The next set of 
tests consider finite-thickness EB beams in more complex geometries. Solutions of these FSI problems are 
obtained using deforming composite grids and the beam solver based on finite elements as discussed in 
Section 7. These problems include the dynamics of a flexible beam separating two fluid chambers, fluid flow 
in a channel with flexible walls, and finally the deformation of a flexible cantilevered beam in a cross fluid 
flow. 


8.1. Verification results for a beam with negligible thickness 

We first consider solutions of two FSI problems involving the interaction of an incompressible fluid and a 
beam of negligible thickness and finite mass. The geometric configuration of both problems is illustrated in 
Figure 6. The fluid domain is given by n(t) = [0,1] x [0, l + rj{x,t)\, where ri{x,t) is the vertical displacement 
of the Euler-Bernoulli (EB) beam which forms the upper surface of the fluid domain. The fluid domain is 
represented by a single transfinite interpolation (TFI) grid which evolves in time according to the computed 
solution for g. The TFI mapping is given by 

{x,y) = g(r,t) = (ri, {I + g{ri,t))r 2 ), r e [0,1] x [0,1]. 

The equations for the fluid are solved in the mapped domain on the unit square with uniform grid spacings 
Ari = Ar 2 . The equations for the EB beam are solved for s = ri S [0,1] using the same grid spacing as the 
fluid grid in the ri direction. 


8.1.1. Manufactured Solutions 

The method of manufactured solutions [32, 33] can be used to construct exact solutions of FSI problems by 
adding forcing functions to the governing equations. The forcing is specified so that a chosen function becomes 
an exact solution to the forced equations. Here the approach is used to verify the stability and accuracy of 
the AMP scheme for the nonlinear FSI problem described above consisting of an initially rectangular fluid 
domain bounded by an EB beam along its top surface. 

The parameters for the fluid are chosen to be p = 1 and fi = 0.05. The thickness of the beam is 
taken to be zero while the finite mass of the beam given by ph is varied to consider the cases of light, 
medium and heavy beams. The remaining parameters of the beam are assumed to he K = 1, T = ph and 
El = Ki = T) =0. The exact displacement of the beam is chosen to be 


ge{x,t) = sin{f^Trx) sm{ftnt), 
T^ft 


(73) 


where a = 0.5 is an amplitude and fx = ft = ‘2 are frequencies. (Note that s = x = ri for this problem.) 
The exact components of velocity and the pressure in the fluid are taken to be 


ui,e(x,t) = -acos{fxTTx)sin{fj;Try)cos{ftTTt), 


^^2,e(x,t) 


a sm{fx;TTx) cos{fxTTy) cos{ftnt) 


dpe 
^ dx 


cos{fj;TTx) sin{fxTry) cos(/t7rt), 


Pe(x, t) = cos{fxTTx) cos{fj;TTy) cos(/t7rt). 


(74) 


where y = y — {1 + pe)- The components of velocity have been chosen so that the continuity equation is 
satisfied in the exact solution, and the vertical velocities of the fluid and the beam match along the interface 
aX y = Pe- Forcing functions are added to the equations for the fluid and structure so that the functions 
in (73) and (74) are exact solutions. 

The stability of the AMP scheme for light beams is illustrated in figure 6 for the case ph = 10“^. The 
image on the left shows the TFI grid for the fluid domain at t = 0.1. The beam lies on the top surface of 
the domain, while the boundary conditions on the other three sides of the domain are taken to be no-slip 
conditions on the left and right sides (v = Vg) and Dirichlet conditions on the bottom (v = and p = Pe)- 
The plot in the center of the figure shows the numerical solution for the vertical component of the fluid 
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velocity, V 2 -, while the plot on the right shows the difference between the computed solution for V 2 and the 
exact solution given by W 2 ,e- Here we observe that the error in V 2 is well behaved in that the magnitude is 
small and it is smooth throughout the domain including the boundaries. The behavior of the error in the 
other components of the solution are similar. 


Grid at £ = 0.1 


V2 at £ = 0.1 


Error in V2 at £ = 0.1 



Figure 6: Manufactured solutions: computational results at t = 0.1 for a light beam {ph — 10 using a grid spacing 
of h — 1/80. Left: deformed grid (coarsened version). Middle: solution component V 2 - Right: error in V 2 - 

A convergence study is shown in Table 7 for light {ph = 10“^), medium {j5h = 1) and heavy {ph = 10^) 
beams using grid spacings Arij = Ar 2 p = hj = l/(10j), j = 1, 2,.... The maximum-norm errors, Ej^'^ for 
solution component q and for grid resolutions j = 1, 2, 4 and 8, are computed for each value of ph. The 
convergence rate is estimated by a least squares fit to the logarithm of the error over the grid resolutions. 
The results in the table for the components of the fluid velocity (^ 1 ,^ 2 ) and pressure p, and for the beam 
displacement ry and velocity pt show the expected second-order accuracy for all three choices of the beam 
mass, and in particular for the case of a light beam when added-mass effects are strong. 

8.1.2. Steady state beam bounding a pressurized fluid chamber 

To illustrate the effectiveness of the AMP coupling scheme in comparison to the traditional FSI coupling, 
we now discuss a problem where the beam deforms as a result of a uniform fluid pressure applied at the 
bottom of the fluid domain. For this problem, an exact solution for the long-time steady state deflection of 
the beam can be found. The steady state can be computed numerically by integrating the time-dependent 
equations from an initially flat beam. It is found that for light beams the traditional scheme is unstable, in 
agreement with the previous stability analysis, while the AMP scheme remains stable and the steady state 
solution is computed accurately. 

For this FSI problem, the material parameters for the fluid are chosen as p = 1 and p = 0.2, while those 
for the beam are chosen as K = 0, T = 5, El = 0, Ki = Ti = 0 and a finite value for ph. The boundary 
conditions on the left and right sides of the fluid domain are given as a no-slip walls (v = 0 ), while a given 
pressure {p = Pa = 5) and zero tangential velocity (ui = 0) are assigned on the bottom of the fluid domain. 
One boundary condition is required at each end of the beam and these are taken as 77 = 0 at each end. (The 
highest spatial derivative in the beam model is two for the chosen parameters.) The initial conditions are 
given as a flat stationary beam (?y = 7 ) = 0) and a fluid at rest. The exact steady-state solution is 

V = 0, P = Pa, P = ^“^2^ 2:e[0,1]. (75) 

As before, ph is varied to illustrate the stability of the time-integration of the equations using the AMP 
coupling scheme for a range of masses of the beam, and also to highlight the instability of the traditional 
FSI coupling scheme for light beams. 

Figure 8 shows the behavior of the displacement and velocity of the beam at its center, a: = 1/2, as a 
function of time for the cases of a light {ph = 10“^), medium {ph = 1) and heavy {ph = 10) beam. For 
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Figure 7: Maximum-norm errors and estimated convergence rates using manufactured solutions. The column labeled 
“r” provides the ratio of the errors at the current grid spacing to that on the next coarser grid. 


each case, there is an early transient stage consisting of oscillations of the beam, but these decay as the 
beam approaches a steady state. Since the damping is supplied by the viscous dissipation in the fluid, the 
oscillations decay more rapidly for lighter beams. At very long times, the solutions for beam displacement 
and velocity approximate the exact steady state solutions. Since the steady state solution is quadratic in 
X, we expect the agreement to be close to exact. The table in Figure 9 gives the max-norm errors in the 
steady state solutions for the light, medium and heavy beams. As expected, the AMP scheme is found to 
be stable and accurate in all three cases. Finally, we note that for this problem the traditional partitioned 
scheme (without sub-time-step iterations) is found to be unstable when ph is less than 5 (approximately). 

8.2. Verification results for a beam with finite thickness 

For the case of FSI problems involving beams with finite thickness, we hrst consider solutions of two 
further benchmark problems to verify the stability and accuracy of the AMP scheme for this more difficult 
configuration. All numerical solutions described in this section, and in the subsequent sections, are computed 
using deforming composite grids and the beam solver based on finite elements. As part of the verification 
tests, we consider both the AMP-PBA and AMP-PBF implementations of the AMP interface conditions. 
For both benchmark problems, the beam has constant nonzero thickness h and separates fluid domains on 
either side. 

8.2.1. Dynamic motion of a flat beam separating two fluid chambers 

We consider a beam separating two initially rectangular fluid domains. The beam is assumed to have 
finite thickness h, and is chosen to initially lie along the a:-axis and occupy the domain f2(0) = [—1,1] x 
[—h/2, h/2]. At time t = 0, the lower fluid domain is the rectangular domain f7_(0) = [—1,1] x [H~, —h/2], 
while the upper fluid domain is n-i-(O) = [—1,1] x [h/2,H'^]. For this first test problem, we assume that 
the beam satisfies a sliding boundary condition ps = EIpsss = 0 at each end, and that the fluid satisfies 
slip-wall conditions along the vertical boundaries. At the top and bottom fluid boundaries, we assume that 
= 0 and p = P±(t), where P±{t) are specified time-dependent pressures. For this configuration, an exact 
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Figure 8: Beam displacement rj and velocity rjt at its center, x = ll2, versus time for light {ph = 10 ®), medium 
{ph = 1) and heavy [ph = 10) beams. 
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Figure 9: Maximum-norm errors at near steady state when t = for light {ph = 10 ®), medium {ph = 1) and heavy 
{ph = 10) beams. 


solution can be found for which the beam moves vertically, but remains flat. In the fluid, the horizontal 
component of the velocity is zero, while the vertical component of the velocity and the pressure depend only 
on y and t. This solution, albeit relatively simple, provides a good test of the stability of the AMP scheme 
for finite-thickness beams and for the two implementations of the AMP interface conditions. 

Working through the equations governing the fluids in each chamber, assuming no dependence on x and 
that the densities in both chambers are equal and given by p, and using the kinematic matching conditions 
at the top and bottom surfaces of the beam, we find that ui = 0 and V 2 = dr]/dt in both fluid chambers, and 


1 

1 rl‘^n 

X S (t) = 

[-1,1] X [H-,r,-h/2], 

t > 0, 


p(x,f) = < 

X G U+(t) = 

[-1,1] X [r, + h/2,H+], 

t > 0. 

(76) 

[ P+(t)-p(p-i7+)_, 


The equation governing an EB beam with all terms involving spatial derivatives set to zero is 

=-KoV-Ki^+p{x,'il-h/2,t)-p{x,ri + h/2,t), a;e[-l,l], t > 0, (77) 

Using the expressions for the fluid pressures in (76), we find that the vertical displacement of the (flat) beam 
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satisfies 


t > 0 , 


(78) 


{-ph + pH)^ + + Kop = AP{t), 

where H = — H~ —h is the total depth of the fluid chambers and AP(t) = P_ (t) — P+ (t) is the difference 

between the applied pressures. We observe that the displacement of the beam satisfies a standard spring 
equation with mass given by the sum of the beam mass, ph, and the fluid added mass, pH, damping given 
by Ki, spring constant given by Kq, and a time-dependent external forcing given by AP{t). Exact solutions 
for r]{t) satisfying (78) for chosen initial conditions are easily found. For example, if 77 = dp/dt = 0 at t = 0, 
and if Ki = 0 and AP{t) is constant, then 

AP 

77 (f) =—{(1 - cos(a;t)}, f> 0 , 
where the frequency of the oscillating beam is given by 


K. 


ph + pH 


Numerical solutions of this ESI problem are computed using the full DCG methodology as discussed in 
Section 7.1. The composite grid ^fp for the geometry is shown in the left plot of Figure 12, and consists 
of a background Cartesian grid and a body-fitted grid adjacent to the beam surface for each fluid chamber. 
The two body-fitted grids are generated independently using the hyperbolic grids generator, and each has 6 
grids lines in the normal direction (so that the grids becomes thinner as the composite grid is refined). The 
grid spacings in each coordinate direction are approximately equal for all component grids and chosen to be 
hj = l/(10j). As before, the index j determines the resolution of the grid. 



Figure 10: Oscillating beam between two fluids for ph — 0.01. Left: beam vertical displacement, 77 , and velocity, 
T]t, compared to the exact solution. Right: errors in 77 and rjt for the AMP-PBF (predict beam internal force) and 
AMP-PBA (predict beam acceleration) schemes. Results are for the coarse grid 

The parameters for the EB beam in the numerical calculations are taken a.sh = 0.1, Kq = 10, f = 10“^, 
El = 10“^, Ki = Tj = 0, and p varied to test light and heavy beams. The parameters for the fluid are 
p = \, p = .02, H~ = —.5 — h/2, and 77+ = .5 -I- h/2 so that H = 1. The pressure difference between the 
top and bottom was taken as AP = 1. Figure 10 shows the computed values of 77 and dp/dt compared to 
the exact solution for a light beam ph = 10“^ on the grid Qfp^- The differences between the computed and 
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Oscillating beam, ph — 0.01 
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Figure 11: Oscillating beam between two fluids. Maximum errors and estimated convergence rates at t = 0.7 
computed using the AMP scheme for a heavy beam, ph = 10 and a light beam ph = 0.01. For this problem the 
convergence rates are close to four since the errors are primarily due to the beam solver which is fourth-order accurate 
in space. 


exact solutions are nearly indistinguishable. Figure 10 also plots the errors in rj and dp/dt for the AMP-PBF 
(predict beam internal force) and AMP-PBA (predict beam acceleration) schemes, and the results are again 
nearly indistinguishable. The tables in Figure 11 shows the max-norm errors and estimated convergence 
rates for a light beam, ph = 10“^, and a heavy beam, ph — t = 0.7. The maximum errors are 

converging close to fourth-order accuracy. Recall that the beam solver is fourth-order accurate in space 
and second-order order accurate in time but converges as fourth-order overall when the time-step scales as 
At (X h^, which is the case for the calculations presented here. 

8.2.2. Steady state beam separating two pressurized fluid chambers 

As a final benchmark test case for which an exact solution is available, we consider the long-time behavior 
of an EB beam separating two initially rectangular pressurized fluid chambers. The geometric configuration 
for this test case is the same as that considered in the previous test for a beam with finite thickness /i, but 
the boundary conditions on the vertical sides are changed so that the beam does not remain flat as time 
evolves. Here we assume no-slip boundary conditions for the fluid on the side walls, and take g = dg/ds = 0 
as clamped boundary conditions for the EB beam. The material properties of the beam are taken to be 
h = 0.1, T = 0, El = 0.2, Ki = 5 and Tf = 0, with p varied to consider light, medium and heavy beams. 
The properties of the fluid are p = 1, p = .05, = 0 and P-{t) = PoR{t), with Pq = 1, where R{t) is a 

smooth ramp function given by 

U35 + i-8A+{70 - 20t)t)t)t\ for 0 < t < 1, 

[l, fort>l. ^ ’ 

The ramp function satisfies R = R' = R” = R'" = 0 at t = 0 and has three continuous derivatives at t = 1. 
Thus, the applied pressure in the lower chamber varies smoothly from P-(0) = 0 to P_(l) = Pq-, and remains 
equal to Pq for t > 1. At t = 0, the fluid in both chambers is at rest, and the displacement and velocity of 
the beam are both zero. 

As in the FSI problem considered in Section 8.1.2, the beam oscillates in response to the nonzero applied 
fluid pressure in the bottom chamber, but then approaches a steady state solution as time increases due to 
the viscous dissipation in the fluid and damping in the beam. At steady state, the velocity in both fluid 
chambers is zero and the fluid pressures are uniform as before, but now for the beam parameters considered 


29 
































here the displacement of the beam is given by 
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2AEJ 
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Figure 12: Composite grids for the partitioned chambers: coarse grid at t = 0 (left) and at t = 10 (right). The 
solution at t = 10 has nearly reached steady state. 

Numerical solutions are computed using the composite grid Q{p\ as before, with grid spacings approxi- 

(“ 2 ) 

mately equal to hj = l/(10j). Figure 12 shows a representative grid, at t = 0 (left) and at t = 10 (right) 
computed for the case ph = 0.1. The solution at t = 10 has nearly reached steady state. The maximum error 
in the beam displacement compared to the steady state solution is approximately 4 x 10“^ at t = 10, and 
reaches 10“^^ hy t = 50. Figure 13 shows the motion of the beam midpoint for different beam densities and 
different algorithmic options with a smaller damping coefficient of Ki = 1. The top-left graph in the figure 
shows the behavior of the displacement, r]{0, t), and velocity, ??t(0, t), as functions of time. Results are shown 
for a light beam, ph = 0 . 01 , and a medium beam, ph = 1, computed on a coarse grid Q^p and a finer grid 

Gfp^- The light beam is excited more rapidly than the medium, and undergoes more rapid oscillations. We 
also observe that the time histories computed using the coarse and fine grid are nearly indistinguishable for 
both beams. The top-right graph in Figure 13 compares results from the AMP scheme and the traditional 
partitioned scheme with sub-iterations (TP-SI) for the case ph = 0.01. Added-mass effects are relatively 
large for this light beam requiring a relaxation parameter in the TP-SI iteration of u = 0.1. The TP-SI 
scheme required an average of approximately 77 sub-iterations per time step to achieve a sub-iteration con¬ 
vergence tolerance of 10“®. The results from the AMP and TP-SI schemes are in good agreement, but the 
computational cost for the TP-SI scheme is significantly larger. 

The dynamic beam problem provides a good test to compare the behavior of the two AMP variations 
described in Section 7.2 for the finite-element based beam solver. It is our experience for this problem, and 
others, that for very light beams the AMP-PBA (predict beam acceleration) algorithm tends to be more 
robust (e.g. requires fewer or no smoothing iterations) than the AMP-PBF (predict beam internal force). 
The AMP-PBF approach requires the evaluation of the third spatial-derivative of the finite element beam 
displacement and this can be sensitive, especially near the ends of the beam where one-sided approximations 
are used (as might be expected). The bottom graph in Figure 13 compares results obtained using the AMP- 
PBA and AMP-PBF for a moderately light beam with ph = 0.1. In this case the results obtained using 
both AMP variations are in excellent agreement. However, for a very light beam, ph = 0.01, the AMP-PBF 
variation has some difficulty and shows poor behavior near the ends of the beam where the acceleration of 
the beam varies very rapidly during the startup phase when the pressure forcing is turned on. It is entirely 
possible that this sub-optimal behavior of the AMP-PBF variation could be rectified with further refinements 
to the algorithm, but we have not pursued this since the alternative AMP-PBA variation behaves well and 
is simpler to implement. 
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AMP-PBA versus AMP-PBF, ph=0.1Q 



Figure 13: Beam under pressure, motion of the beam mid-point. Top left: r) and rjt for a light {ph — 0.01) and 
medium {ph — 1) beam using a coarse grid (colored curves) and a finer grid (black curves). Top right: AMP 
versus the TP-SI scheme for the light beam, ph — 0.01. Bottom: A comparsion of the two AMP schemes for coupling 
to the finite-element beam solver on grid for ph = 0.1. 


8.3. Traveling-wave pressure-pulse in an elastic tube 

The propagation of a pressure pulse through a two-dimensional elastic tube can be used as a model for the 
flow of blood in a large artery or vein following the problem proposed in [34], and later considered in [13, 21] 
to study the stability of partitioned schemes. Added-mass effects are important in this application since 
the density of the elastic wall is close to that of blood and also since the tube is long and thin. Traditional 
partitioned methods require many sub-iterations per time step to remain stable for this problem. The AMP 
algorithm is shown to provide stable results with no sub-iterations per time step. 

The domain for the fluid is the channel n(t) = [0, L] x [0, R{t)], where L is the length of the channel and 
R{t) is the distance from the axis of symmetry of the channel at y = 0 and the elastic wall. We consider the 
case when L — Q and i?(0) = 0.5. The boundary conditions for the fluid domain are as follows. The bottom 
boundary is a slip-wall where V 2 = 0 and dvi/dy = 0. The top boundary is the flexible beam where the 
AMP interface conditions are applied. On the left boundary the tangential component of the velocity is set 
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Figure 14: Composite grids for the traveling-wave pressure-pulse in a flexible tube: coarse grid at times t — 1.0, 
1.5 and 2.0. 


to zero, V 2 = 0, and the pressure is specified to be the time-dependent pulse 


p{0,y,t) 


Pmax Sin(7rt/ti,;iax) for 0 ^ ^ tmax; 

0 for t > tinax; 


(80) 


where Pmax and tmax determine the magnitude and duration of the pulse (values given below). The right 
boundary is an outflow boundary where p = 0 and the components of velocity are extrapolated. 

The composite grid for the fluid domain at various times is shown in Figure 14. The composite grid 
consists of a background Cartesian grid (shown in blue) together with a hyperbolic grid adjacent to the top 
interface (shown in green). The hyperbolic grid is generated with 6 grids lines in the normal direction and 
thus this boundary-fitted component grid becomes thinner as the composite grid is refined. The grid spacing 
is chosen to be hj = 1/(10_)) approximately. 

We perform simulations of the elastic tube problem for two choices of the parameters of the fluid and 
the beam, identified as Parameter Sets I and II. The parameter sets are based on the parameters chosen in 
Fernandez et al. [21], but with some differences as described below. For both choices, the beam is taken 
to be a generalized string model with El = 0 (i.e. no fourth-order spatial derivative term in the beam 
equation). The AMP-PBA (predict-beam-acceleration) variation is used for all calculations in this section 
and the projected interface velocity is smoothed with three iterations of the fourth-order filter, as described 
in Section 7.3. 


8.3.1. Parameter Set I 

Dimensionless parameters for this set are taken to be p = 1, /r = 10“^, Pmax = 2 and tmax = 0.75 for the 
fluid, and ph = 0.11, T = 5/6 « .83333, Kq = 40/3 « 13.333, Ki = ph/100 and Ti = 5/60 « .08333 for the 
beam. The term with coefficient Kq acts as a linear restoring force while the term with coefficient Ti is a 
visco-elastic damping term that tends to smooth high-frequency oscillations in space. The parameters used 
for this first set are slightly different than the ones used in [21]. Our parameters give a larger amplitude 
of the deformation of the beam and a thicker boundary layer in the fluid near the surface of the beam. 
A larger deformation of the beam is a more severe test of our DCG approach, and a thicker boundary 
layer enables grid-convergence tests demonstrating second-order accuracy using coarser grids. (The choice 
of viscosity used in [21] requires a very fine grid to resolve the very thin boundary layer, see Parameter Set 
II in Section 8.3.2 below). We also note that the fluid tractions and velocities are transferred directly on the 
beam reference curve for the simulations described in [21], whereas our approach applies these quantities on 
the beam surface as described in Section 7.2. In order to approximate a beam of negligible thickness with 
our approach, we choose a small value for the beam thickness ofh = 10“^. 
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Figure 15: Parameter set I: Traveling-wave pressure-pulse in a flexible tube. Contours of p, vi and U 2 at t = 1.5 

/•q’\ 

computed on grid For clarity, the top surface curve is also is shown in red with the displacement amplified by 

a factor of 5. 


Figure 15 shows shaded contours of the fluid pressure and the components of the fluid velocity at t = 1.5. 
The pressure pulse at this time has created a bulge in the flexible channel near its half-length, and for clarity 
the centerline-curve of the beam on the top surface is drawn as a red curve with the displacement amplified 
by a factor of 5. The traveling pressure pulse is seen to consist of a localized high pressure region where 
the beam displacement is positive. The pressure is largest at the surface near the maximum displacement of 
the beam. The horizontal component of the velocity is generally large and positive within the pulse with a 
boundary layer near the top surface to match the no-slip condition there. The vertical component of velocity 
is positive ahead of the pulse and negative behind it confirming that the pulse is traveling in the positive xi 
direction. 

The top left plot in Figure 16 shows the evolution of the shape of the surface at four different times. 
The surface is initially displaced upwards on the left side due to the high pressure at inflow (80) which 
acts over the time interval [0,tniax]; where tniax = 0.75. The surface pulse then propagates to the right and 
slowly decays in amplitude. The top right plot shows the beam displacement at t = 1.5 computed using 
grids with increasing resolution. The solutions are seen to converge with the curves on two finest grids being 
nearly indistinguishable. The two plots on the bottom of Figure 16 show the grid convergence of the vertical 
component of the fluid velocity and the fluid pressure at the surface of the beam at t = 1.5. Here, we also 
observe excellent convergence of these two fluid variables. 

To provide a more quantitative estimate of the solution accuracy, a self-convergence study is performed. 
Given a sequence of three or more grids of increasing resolution, a posteriori estimates of the errors and 
convergence rates can be computed using the Richardson extrapolation procedure described in [35]. These 
self-convergence estimates assume that the numerical solution converges to a limiting solution, but does 
not assume that some chosen fine-grid solution is the exact solution. A posteriori estimates computed in 
this way are given in the table in Figure 17 for four grids of increasing resolution (the convergence rate 
computation uses the three finest resolutions). The table provides estimates of the max-norm errors for 
various components of the solutions and the corresponding estimates of the convergence rates. We observe 
that all solution components are converging with approximately second-order accuracy. Note that due to the 
manner in which rates are computed, the ratios of errors in the columns headed “r” are the average ratios 
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Flexible channel, surface displacement 



Flexible channel, surface displacement 




Figure 16: Traveling-wave pressure-pulse in a flexible tube. Top left: surface displacement over time showing the 
formation and propogation of the pulse. Top right: surface displacement at t = 1.5 for grids k = 2,4,8,16 
showing the convergence of the surface as the grid is refined. Bottom left and lower right: convergence of V 2 and p, 
respectively, as the grid is refined. 


and thus the same for each grid resolution. 


Traveling wave pressure-pulse 
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Figure 17: Traveling wave pressure-pulse (parameter set I): Maximum errors and estimated convergence rates at 
t = 1.5 computed using the AMP scheme. 

Figure 18 compares results from the AMP scheme with the traditional-partitioned scheme with sub-time- 

/n\ 

step iterations (TP-SI) on the coarse grid Get - The number of sub-iterations per time step required for the 
TP-SI scheme is approximately 15 on average, although near the start of the run almost 250 iterations are 
needed. (A convergence tolerance of 10“^ on the sub-iterations is used.) A simple under-relaxed fixed-point 


34 









































Flexible channel, surface displacement 



Figure 18: Traveling-wave pressure-pulse in a flexible tube. A comparison between the AMP scheme and the TP 

( 2 ) 

scheme with sub-time-step iterations at t = 1.5 with grid ■ The surface displacement, velocity V 2 and pressure are 
shown (with markers placed every 8 grid points). The results of the two schemes are seen to be nearly indistinguishable. 


iteration is used with the value for the under-relaxation parameter taken as a; = 0.025 (the fact that such a 
small value is needed indicates that the problem has large added-mass effects). The solutions from the AMP 
and TP-SI are in excellent agreement. 

8.3.2. Parameter Set II 

The parameter values for this set match those of Fernandez et al. [21] so that a direct comparison of the 
results can be made. We prefer to work with dimensionless parameters obtained by scaling with the reference 
length Lq = 1 cm, velocity Uq = ItP cm/s, density po = 1 gm/cm^, and time Tg = Lq/Uq = 10“^ s. These 
scalings give the corresponding dimensionless parameters p = 1, p = 3.5 x 10“"^, Pmax = 2 and tmax = 0.5 
for the fluid, and ph = 0.11, T = 2.5, Kq = 40, Ki = 1.1 x 10“^ and Ti = 0.25 for the beam. As in the 
previous set of parameters, the product ph = 0.11 is fixed but h = 10 “^ and p = 110 are chosen to reduce 
the effects of a finite-width beam. 

Figure 19 compares the surface displacement of the beam at t = 1.5 computed with different grids using 
the AMP-PBA algorithm. The composite grids used in these calculations are the same as those used in the 
previous section. We observed that the surface displacement appears to be converging in a manner consistent 
with second-order accuracy. The figure also shows the results from Fernandez et al. [21] using a first-order 
accurate monolithic scheme on a grid with spacing h = 1/320 (equal to the spacing on grid '). Although 
the monolithic results do not exactly match the fine grid AMP results, we note that the monolithic results 
do not yet appear to be grid converged since they differ by a fair amount from the monolithic results on a 
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Figure 19: Paramater set II: Traveling-wave pressure-pulse in a flexible tube. Surface displacement as the grid are 
refined. Results for monolithic scheme are from Fernandez et.al. [21]. 


grid coarsened by a factor of 2, see [21]. 



Figure 20: Parameter set II: Traveling-wave pressure-pulse in a flexible tube. Contours of p, v\ and U 2 at t = 1.0 

/o') 

computed on grid ^ For clarity, the top surface curve is also is shown in red with the displacement amplified by 
a factor of 10. 


Figure 20 shows shaded contours and surface displacements of the solution at t = 1.0 computed using 
the AMP algorithm. We note that the behavior of the solution is qualitatively similar to that observed 
previously for Parameter Set I. The contour plot of the horizontal component of velocity, Ui, however, shows 
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the very thin boundary layer that forms near the interface. 


8.4- Beam in a cross flow 

As a final demonstration of our AMP algorithm, we consider flow past a flexible beam that partially 
blocks a fluid channel as shown in Figure 21. A flow from left to right causes the beam to initially bend to 
the right and then oscillate. Eventually the flow approaches a steady-state consisting of a long recirculation 
region downstream of the deflected beam. The computed solution from the AMP scheme for a light beam 
is compared with the solution obtained using the TP-SI scheme. Results for beams of different densities are 
also shown and used to estimate the added mass for this configuration. 

The geometry of the problem is illustrated in Figure 21 and consists of an initially vertical beam, 
with length 1 = 1, that extends into a rectangular fluid channel of height yt, = 2 and overall dimensions 
[Xa,Xb] X [uajUb] = [”3,8] X [0,2]. The thickness of the beam at its clamped base is ^(0) = 0.2 (the beam 
reference line is located along the centerline) and this remains approximately constant along the beam until 
its rounded free end^. Note that while the beam reference-line (red curve in Figure 21) has an undeformed 
length of Z = 1 , the actual distance from the base of the beam to the rounded tip at the top is slightly longer®. 
The base of the beam is at the origin of the computational domain, and the upstream vertical boundary 
of the channel is at Xa = —3 while the downstream boundary is at Xb = 8 . The grid for this geometry, 
denoted by Gflcfl shown previously in Figure 5, and consists of a background Cartesian grid in the fluid 
channel together with a boundary-fitted hyperbolic grid that extends around the beam. The hyperbolic grid 
is generated with 6 grids lines in the normal direction, and thus this grid becomes thinner as the composite 
grid is refined as in previous examples. The grid spacing is chosen to be approximately hj = l/(10j). 

The parameters for the fluid are p = 1 and p = 0.05, while the parameters for the beam are El = 5, 
K = T = Ki = Ti = 0, and with a range of values for the beam density p as noted below. The conditions 
on the boundaries of the fluid channel are chosen as a no-slip wall on the bottom boundary at y = 0, a 
slip-wall on the top boundary at y = y^, inflow on the left at a: = a;a, and outflow on the right at x = Xb- 
The outflow boundary condition is a Robin condition on p and extrapolation of the velocity components. 
The inflow velocity profile is taken to be uniform over most of the inflow boundary but parabolic near the 
bottom to match the no-slip boundary condition on the lower wall. The inflow profile is given as 


vi{xb,y,t) = Vifly)R{t), Vifly) 


kmax for y ^ d, 

Kiax(y/d)^ for y <d, 


(81) 


where Vmax = 1, d = 0.2 and R{t) is the ramp function given previously in (79) so that the horizontal 
component of the inflow velocity transitions smoothly in time from = 0 at t = 0 to f i = Idn(y) for t > 1. 
The vertical component of the inflow velocity at x = Xa is zero. 

Figure 22 shows the behavior of the displacement and velocity of the tip of the beam for the case of a 
light beam with p = 1. Solutions are computed using the AMP scheme with composite grids j = 2,4, 8 , 
to indicate the grid convergence. We also show the behavior of the displacement and velocity computed 
using the TP-SI scheme for a fine grid Gbc}- zoomed views show excellent grid convergence for the 
AMP scheme, and excellent agreement with the solution obtained using the TP-SI scheme. We note that 
the TP-SI scheme required, on average, approximately 19 sub-iterations per time-step using grid and 

6 sub-iterations using grid Gb^'f, with an iteration tolerance of 10 “® in both cases. 

Figure 23 shows the tip displacement and velocity for beams of different densities computed using the 
composite grid Gbt'f- We observe damped oscillations for all beams, with the damping for the heavy beam 


^ The reason that the beam at the tip is rounded off, rather than having sharp corners, is so that a high-quality 
grid can be constructed that resolves the rapid variation of the solution around the tip. In addition, the solution 
and its derivatives will remain smooth near the rounded tip and this will give greater confidence that the numerical 
solution is converging to the true solution in this region. 

® The actual boundary and rounded end of the beam, and hence its thickness, is defined in terms of hyperbolic 
trigonometric fnnctions through the SmoothedPolygon mapping [36]. The precise definition of the grid along with 
scripts to run the simulations presented here will be available with the Overture software at overtureFraniework.org. 
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Figure 21: Light beam with p = 1 in a cross flow: Streamlines of the flow near steady state at t = 20 computed on 
the composite grid 



Figure 22: Light beam in a cross flow: tip displacement and velocity as the grid is refined for parameters El = 5, 
p= 1. 


being the weakest and its approximate period (time between successive peaks) being the largest. We also 
note higher-frequency sub-harmonic motion of the beam tip, especially for the heavy beam. All solutions 
eventually approach a steady state where the displacement of the beam tip is ptip ~ 0.075 independent of 
the beam density. 

In the absence of the fluid, the frequency, wq, and period, tq, of the lowest mode of a cantilevered beam 
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Beam in a cross flow: tip displacement Beam in a cross flow: tip velocity 



Figure 23: Beam in a cross flow: comparison of beam tip displacement and velocity for beams of different densities. 
The solution was computed on grid with El = 5. 


are well known [37] and given by 


Wo 


( 1 . 875)2 




(82) 


By comparing the computed frequency and period of the damped oscillating beam from the FSI simulation, 
the added-mass of the fluid can be estimated. In view of the analysis discussed in Section 6 we expect that 
the frequency of the beam, when accounting for added mass effects, would be approximately given by 




(1.875)2 


El 


{ph + Ma)l'^ 


Tam 


27r 

WAM ’ 


(83) 


where Ma is the added mass of the fluid (see equation (47)). Thus, using (82) and (83), we have the estimate 


Ma = ph{{TAM/To)‘^ - l}. (84) 

For p = 100, h = 0.2, I = 1 and El = 5, we find tq « 3.57 which is close to the period estimated from 
Figure 23 of tam ~ 3.6. Using (84), we find Ma/{ph) « .02, which is small as expected for this heavy beam. 

For a lighter beam with p = 1, we obtain tq « 0.357 and tam ~ -99, which gives Ma « 1.3, and this 
value is approximately 6.5 times larger than the beam mass of ph = 0.2. Based on the earlier mode analysis 
and the stability result in (63), we can estimate that the traditional scheme would become unconditionally 
unstable for p = Ma/h = 1.3/0.2 = 6.5 for At approaching zero. Computations of the full problem indicate 
that the TP scheme actually becomes unstable at a larger value of p « 13. This larger value is expected 
since the time-step for this computation is chosen based on the stability of the domain solvers, as might be 
done in practice, and not reduced to account for added-mass effects. 


9. Conclusions 

This article makes a significant contribution to resolving an important and long-standing difficulty for 
fluid-structure interaction (FSI) problems involving incompressible fluids and thin structures. Such FSI 
problems arise in many applications in science and engineering, and their accurate solution is an important 
and active field of research. In a pair of recent papers [4, 5] , we have developed new added-mass partitioned 
(AMP) algorithms for linearized FSI problems coupling incompressible (Stokes) flows and linearly elastic bulk 
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solids and thin structures. An important feature of the AMP algorithms is that they are stable, without sub- 
time-step iterations, independent of added-mass effects. The present paper represents a first and important 
extension of the AMP algorithm in [5] to fully nonlinear FSI problems for the case of thin structures. The 
AMP interface conditions and time-stepping algorithm has been fully described for this class of FSI problems. 
The key AMP interface condition is a Robin (mixed) boundary condition for the fluid pressure, and for a 
general beam this becomes a generalized Robin condition that couples points on the circumference of the 
beam cross-section. For a two-dimensional beam with fluid on two sides, for example, this condition involves 
the pressure on opposite sides of the beam. While the AMP algorithm has been described for general beams, 
the focus has been on beams in two dimensions described by the Euler-Bernoulli (EB) equations. A normal¬ 
mode stability analysis has been performed for a linearized problem involving an EB beam separating fluid 
domains on either side. The analysis extends that given in [5] and shows that the AMP scheme is insensitive 
to added-mass effects for this more general problem configuration, while a traditional partitioned scheme 
(without sub-time-step iterations) becomes unconditionally unstable when added-mass effects are sufficiently 
strong. 

The general numerical framework for handling finite-amplitude deformations of the structure is based 
on the use of moving and deforming composite grids, previously used in the context FSI problems involving 
compressible flows coupled to rigid solids [30, 38], and also to linear [6] and nonlinear [7] elastic solids. The 
equations governing the incompressible fluid are solved using a fractional-step finite-difference scheme on 
deforming grids, while the equations for the EB beam are solved using either finite differences or a finite- 
element approach. For the case of the finite-element beam solver, special treatment of the AMP interface 
conditions are required in order to couple to the finite-difference fluid solver and we have described two 
approaches for implementing this coupling. 

An important feature of this paper is the use of benchmark problems of increasing complexity to illustrate 
the behavior of the new AMP algorithm. In addition to carefully verifying the stability and second-order 
accuracy of the present algorithm, the results of the benchmark problems can be useful to other researchers 
developing other FSI algorithms. The AMP scheme was first verified using manufactured solutions for the 
case of a deforming beam with a fluid domain on one side. The solution was next computed for a problem 
(with known steady state solution) of a deformed beam next to one or two pressurized fluid chambers, and 
for a time-dependent oscillating flat beam (with known solution). Results were obtained using a traditional- 
partitioned (TP) scheme with sub-time-step iterations and these agreed with those obtained using the AMP 
scheme. The AMP scheme was also applied to a problem from the literature for a simplified model of 
an artery: solutions were computed for a traveling wave pressure-pulse propagating through a deforming 
channel. Grid convergence studies were performed and a comparison was made to prior results available in 
the literature. Lastly the motion of a cantilevered beam in a fluid channel was considered, and solutions for 
beams with different densities were compared and these results were used to estimate the added mass. In all 
benchmark cases considered the numerical solutions from the AMP algorithm remained accurate and stable 
even for the case of light beams when added-mass effects are large. 

In future work we will consider more general nonlinear beam models as well as FSI problems involving 
beams, tubes, and shells in three space dimensions. 
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